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Using an unbiased quantum Monte Carlo method, we obtain convincing evidence of the existence 
of a checkerboard supersolid at a commensurate filling factor 1/2 (commensurate supersolid) in the 
soft-core Bose-Hubbard model with nearest-neighbor repulsions on a cubic lattice. In conventional 
cases, supersolids are realized at incommensurate filling factors by a doped-defect-condensation 
mechanism, where particles (holes) doped into a perfect crystal act as interstitials (vacancies) and 
delocalize in the crystal order. However, in the above model, a supersolid state is stabilized even at 
the commensurate filling factor 1/2 without doping. By performing grand canonical simulations, we 
obtain a ground-state phase diagram that suggests the existence of a supersolid at a commensurate 
filling. To obtain direct evidence of the commensurate supersolid, we next perform simulations in 
canonical ensembles at a particle density p = 1/2 and exclude the possibility of phase separation. 
From the obtained snapshots, we discuss the microscopic structure and observe that interstitial- 
vacancy pairs are unbound in the crystal order. 

PACS numbers: 03.75.Hh, 05.30. Jp, 67.85.-d 



Since the observation of a nonclassical moment of in- 
ertia in solid 4 He by Kim and Chan[l], Q, the possible 
existence of supersolid has been discussed intensively. 
Regarding the possibility of a bulk supersolid in 4 He, 
Prokof'ev and Svistunov pointed out that, in the ab- 
sence of symmetry between vacancies and interstitials, 
a commensurate supersolid in continuous space has zero 
probability of being observed in natureQ. Thus, they 
suggested that another scenario should be considered to 
interpret the decrease in nonclassical rotational inertia 
observed in 4 He. In contrast to the case of continuous 
spaces, the presence of supersolid states in lattice systems 
has been established by unbiased quantum Monte Carlo 
simulations 043, particularly in the case when particles 
are added to or removed from perfect crystals. As a re- 
sult of recent progress in experiments on cold atoms and 
molecules, optical lattice systems have become highly 
promising systems for realizing the supersolid state. In 
conventional cases, we can call such a supersolid state an 
incommensurate supersolid state, because it is realized 
by doping defects such as interstitials or vacancies into 
perfect crystals. 

Although a commensurate supersolid in continuous 
space has zero probability of being found in nature, this 
does not apply to systems with explicitly broken trans- 
lational symmetry such as lattice systems Q. Among the 
optical lattice systems, one of the most promising can- 
didates for realizing a supersolid is a checkerboard-type 
supersolid near filling factor 1/2, because it can be real- 
ized only in the presence of appropriate nearest-neighbor 
repulsions of soft-core bosons on simple square or cubic 
lattices @, El 03 • Therefore, it is an interesting ques- 
tion whether such a checkerboard supersolid can be sta- 
bilized even at the commensurate filling factor 1/2 in 
the absence of doping. Regarding the soft-core extended 
Bose-Hubbard model with nearest-neighbor repulsions, 



the mean-field analysis [TlJ and the Gutz wilier variational 
method [HI, [r| predicted that a checkerboard supersolid 
is stabilized at the commensurate filling factor (particle 
density) p = 1/2 in addition to above and below p = 1/2. 
In contrast to these studies, exact quantum Monte Carlo 
studies have found evidence of a supersolid phase only 
above p = 1/2 in ID and 2D systems @, @|. However, 
recent quantum Monte Carlo studies have revealed that 
a supersolid phase exists below p = 1/2 as well as above 
p = 1/2 in the case of a 3D cubic lattice [TJ E3- This is 
consistent with the prediction by the mean-field analysis 
and the Gutzwiller variational method, and also suggests 
the presence of a supersolid phase at p = 1/2. To obtain 
direct evidence of a supersolid existing at the commen- 
surate density p = 1/2, it is necessary to consider the 
canonical ensemble at p = 1/2 and exclude the possi- 
bility of a phase separation, neither of which have been 
investigated previously. In this Letter, using an unbi- 
ased quantum Monte Carlo method, we show convinc- 
ing evidence of a commensurate supersolid existing in 
the above-mentioned 3D system. The finite-temperature 
transition to a commensurate supersolid is also investi- 
gated. Finally, we show snapshots of the commensurate 
supersolid and discuss a microscopic structure. 

The model considered here is the soft-core Bosc Hub- 
bard model with nearest neighbor repulsions on cubic 
lattices. The Hamiltonian is given by 
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Here, b\(bi) and m are the bosonic creation (annihilation) 
operator on site i and the particle number operator de- 
fined as rii — b\bi, respectively. Furthermore, t, p, U, and 
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FIG. 1: (Color online) Ground-state phase diagram of the 
soft-core Bose Hubbard model with nearest neighbor repul- 
sions on cubic lattices. Error bars are drawn but most of 
them are much smaller than the symbol size (here and in the 
following figures). Lines are used to guide the eyes. The black 
points at t/U = and the black solid line separating p = 
from SF are determined analytically. 



V represent the hopping parameter, the chemical poten- 
tial, the on-site repulsion, and the nearest-neighbor re- 
pulsion, respectively. The summation is taken over 
all pairs of nearest-neighbor sites. In our simulations, we 
treat N = L 3, systems with the periodic boundary condi- 
tion. To investigate the above model, we performed un- 
biased quantum Monte Carlo simulations based on a hy- 
brid algorithm of a worm algorithm (l~7l - [l9| and an O(N) 
Monte Carlo method[20l|. This hybrid algorithm enables 
efficient simulation even in the presence of off-site inter- 
actions between bosons. 

We show the ground-state phase diagram of the model 
in Fig. [TJ In our simulations, we set the nearest-neighbor 
repulsion as V/U = 1/z, where z = 6 is the coordination 
number. For < p < 1, we confirmed the presence of a 
supcrfluid (SF) phase, a checkerboard (CB) solid phase 
at p = 1/2, a Mott insulator (MI) phase at p = 1, and 
a supcrsolid (SS) phase. The phase diagram is quali- 
tatively in agreement with previous results obtained by 
the mean- field analysis [ll| and the Gutzwiller variational 
method(l5l.[2l|. Since the added particles on the checker- 
board solid are subjected to an almost flat potential on 
the checkerboard-type background at U/V ~ 1/z, they 
can gain high kinetic ener gy a nd delocalize without de- 
stroying the crystal order 11| . Consequently, a broad 
interstitial-based SS phase appears above the CB solid 
phase. 

To determine the phase boundaries, we calculated the 
particle density p = l/N^Arii), the supcrfluid stiff- 
ness p, = (W 2 )T/(Lt), and the structure factor S(k) = 
^lNY J ,,e lk - r ^{(n l n j )~{n l ) 2 ). Here, W, k, and r tj are 
the winding number vector in the path integral represen- 
tation, the wave vector, and the relative position vector 
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FIG. 2: (Color online) (a) Physical quantities as functions 
of the chemical potential at t/U = 0.05. (b) Determination 
of the quantum critical point between SS and SF from the 
intersection of the Binder ratios for different system sizes, 
(c) Physical quantities as functions of the chemical potential 
at t/U = 0.057. (d) Physical quantities as functions of the 
particle density for the same parameter set as that in (c). 



between sites i and j, respectively. Furthermore, (■ • ■ ) 
indicates the thermal expectation value. In Fig. HJa), 
the results are shown as functions of the chemical po- 
tential p/U at (t/U,V/U,T/t) = (0.05, 1/z, 0.3). Since 
the broken symmetries in CB and SF are not associated 
with each other, it is expected that a first-order phase 
transition takes place. We estimated its phase bound- 
ary from the position of the discontinuity of p, p s and 
S (it , tt , ir) . As a simple but reasonable way to deter- 
mine the boundary between CB and SS, we calculated 
the zero-momentum Green function and estimated the 
energy gap for inducing the particle/hole excitation{22l]. 
In the same manner, the boundary of MI is also obtained 
from the zero-momentum Green function. The determi- 
nation of the phase boundary between SS and SF requires 
a more qualitative consideration. As shown in Fig. EJb), 
wc estimated it from the intersection of the Binder ratio 
g = l/2[3— (m 4 )/(m 2 ) 2 ] for different system sizes, where 

m = l/NJ2i n i el ^ Vi at fc = (7T, 7r,7r). In this analysis, 
we assumed that the dynamical critical exponent is equal 
to 1, and we fixed the temperature as T/t oc 1/L. For 
t/U < 0.03, we observed a discontinuity in the quanti- 
ties at the boundary between SS and SF. Therefore, we 
determined the boundary from the position of the dis- 
continuity. 

For t/U > 0.055, the SS phase covers the tip of the 
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FIG. 3: (Color online) (a) Finite-temperature transition to 
the commensurate supersolid state, (b) Finite-size scaling 
of S(tt, 7r, it) for the normal-to-CB transition, (c) Finite-size 
scaling of p a for the CB-to-SS transition. 



CB lobe. This result suggests the possible existence of 
commensurate supcrsolids, because we can expect that 
the system with a fixed density p = 1/2 undergoes a 
phase transition from the CB solid state not directly to 
the SF state, but to the SS state as t increases. This 
can be also expected from p dependence of the mea- 
sured quantities. In Fig. HJc), the physical quantities 
at (t/U,V/U,T/t) = (0.057, l/z,0.3) are shown as func- 
tions of the chemical potential. As the chemical potential 
increases, the SF-to-SS transition occurs at p < 1/2, and 
the particle density increases to above p — 1/2. Since 
there is no insulating phase in which the particle density 
has a plateau and the superfluid stiffness vanishes, we 
expect that the SS phase remains even at p = 1/2. To 
explain this in more detail, we also show p s and >!5(7r, 7r, n) 
as functions of the particle density p in Fig. [5{d). These 
results are obtained by making bins for particle densities 
with a finite width and categorizing each sample into 
the appropriate bin. Although these results support the 
hypothesis that both finite p s and S(ir, 7r,7r) survive at 
p = 1/2, we also have to exclude the possibility of a phase 
separation at exactly p = 1/2, to conclude the presence 
of a supersolid at the commensurate filling factor 1/2. 

To obtain direct evidence of the existence of a com- 
mensurate supersolid, we next perform simulations for 
p = 1/2. The following results were obtained by using 
only samples with p that is exactly equal to 1/2. In other 
words, we performed quantum Monte Carlo calculations 
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FIG. 4: (Color online) Finite-temperature phase diagram at 
the fixed density p — 1/2. Circles and squares denote the criti- 
cal temperatures for the checkerboard order and superfluidity, 
respectively. The disordered phase at finite temperatures is 
referred to as normal liquid (NL) phase. 



in a canonical ensemble at p = 1/2. In Fig. [31 the tem- 
perature dependences of p s and S(7r, 7r, 7r) are shown at 
(t/U,V/U) = (0.057, 1/z). We can confirm finite values 
of p s and S(n, 7r, 7r) at low temperatures. In our sim- 
ulations, p s and S(ir, 7r, 7r) for all replicas, which were 
simulated under different initial conditions and with dif- 
ferent seeds of random number, converged to the same 
values. This indicates the absence of a phase separation. 

To clarify the existence of a finite-temperature transi- 
tion to the commensurate supersolid state, we performed 
finite-size scaling analysis. Since the breaking symme- 
tries in the supersolid state are and U(l), which are 
related to the checkerboard order and the superfluidity, 
respectively, two successive transitions of the Ising-type 
universality class and the XY universality class are ex- 
pected. We analyze the finite-size-scaling behavior of the 
structure factor and the superfluid stiffness by consid- 
ering the scaling forms S(Tr,ir,n)L 1+ri = /(SL 1 ^) and 
p s L = g{5L 1 / u ), respectively Here, S is defined as 
5 = (T — T c )/T c , and T c is the critical temperature. Us- 
ing the critical exponents of the 3D Ising universality 
class (i/=0.63002 and ?7=O.O36270) for the structure 
factor and the 3D XY universality class (V=0.67155[24j|) 
for the superfluid stiffness, we successfully observed data 
collapses. As shown in Figs. EJb) and (c), the data of 
system sizes larger than L ~ 16 are necessary to obtain 
good scaling results. These results of finite-size scaling 
also strongly indicate the absence of a phase separation 
at p = 1/2. Therefore, we conclude that a true com- 
mensurate checkerboard supersolid exists in the present 
model. In Fig. 01 we summarize the transition points in 
the finite-temperature phase diagram at p = 1/2. 

Finally, we discuss the microscopic structure of the 
commensurate supersolid. As discussed in Ref. Q, in 
the supsersolid state, interstitials and vacancies do not 
form bound pairs and can be found arbitrarily far from 
each other. The delocalized defects give rise to super- 
fluidity in crystal orders. In contrast, in the insulating 
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FIG. 5: (Color online) Cross-sectional snapshots 
at a fixed z-coordinate and imaginary time in 
N—20 3 systems. (a) Checkerboard solid state at 
(t/U,V/U,T/U)=(0.052, 1/2,0.02) and (b) commensu- 
rate supersolid state at (t/U, V/U, T/U)=(0.056, 1/ z, 0.02). 
Particle filling is exactly p = 1/2 in both cases, although this 
cannot be confirmed from a single cross section above. Sites 
are denoted as rectangles. Empty rectangles, blue rectangles, 
and red rectangles indicate empty sites, singly occupied sites, 
and doubly occupied sites, respectively. In (a), the green 
rectangular frame represents a bounded interstitial-vacancy 
(particle-hole) pair confined within a nearest-neighbor site. 
In (b), the existence of several doubly occupied sites suggests 
the presence of interstitials and vacancies that are unbound 
beyond the nearest-neighbor sites. 



phases, locally created particle-hole pairs are typically 
confined within nearest-neighbor sites, as has been well 
discussed for the Mott insulator at p = 1 [25[ . To qualita- 
tively confirm the characteristic features, we show snap- 
shots at the commensurate density p = 1/2 in Fig. [5] 
In the checkerboard insulating phase [Fig. EI a)], we ob- 
serve bound interstitial-vacancy pairs that are typically 



confined within nearest-neighbor sites. In contrast, in 
the commensurate supersolid state [Fig. E{b)], we can 
easily observe interstitials and vacancies that are sepa- 
rated by a distance exceeding the nearest-neighbor dis- 
tance. This can be confirmed from the presence of dou- 
bly occupied sites because such sites are never created 
from interstitial-vacancy pairs confined within nearest- 
neighbor sites. Since we did not dope any interstitials or 
vacancies into the commensurate solid, the defects orig- 
inated from the unbound interstitial-vacancy pairs. In 
this sense, it is expected that the mechanism of the CB- 
to-SS transition is similar to the released doublon-holon 
mechanism of the MI-to-SF transition at p = 1 25| . 



In conclusion, using exact quantum Monte Carlo sim- 
ulations, we obtained direct evidence of the existence fo 
a checkerboard supersolid at p = 1/2 in 3D bosonic sys- 
tems. By obtaining a finite-temperature phase diagram 
at p = 1/2, we also showed the range of parameters in 
which this supersolid exists, including its temperature. 
Finally, from snapshots, we qualitatively discussed the 
its microscopic structure. Although the snapshots sug- 
gest the existence of unbound interstitial-vacancy pairs 
in the commensurate supersolid, further investigation is 
needed to confirm this quantitatively. 
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